Languguage OS 2/Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO/language/awe/awe-full.lha/Awe2/DoNotUseThisSrc/NewRandom

Text File | 1989-11-29 | 47KB | 1,588 lines

(Message inbox:3183)
Return-Path: [e-mail]
Received: from truth.waikato.ac.nz by june.cs.washington.edu (5.61/7.0j)
    id AA06839; Thu, 31 Aug 89 19:38:53 -0700
Message-Id: <[e-mail]>
Received: from wnv.dsir.govt.nz by waikato.ac.nz; Fri, 1 Sep 89 14:07 +1200
Date: Fri, 1 Sep 89 14:09 NZST
From: ROBERT <[e-mail]>
Subject: RE: Re: random number package
To: [e-mail]
X-Vms-To: IN%"[e-mail]",SRWMRBD

RANDOM NUMBER GENERATOR PACKAGE


Copyright (C) 1989: R B Davies and DSIR

Permission is granted to use or distribute but not to sell.

The gamma function routine is adapted from "Numerical Recipies in C" by Press,
Flannery, Teukolsky, Vetterling published by the Cambridge University Press.

This is a package for generating sequences of random numbers from a wide
variety of distributions. It is particularly appropriate for the situation
where one requires sequences of identically distributed random numbers since
the set up time for each type of distribution is relatively long but it is
very efficient for generating each new random number. The package includes
"classes" for generating random numbers from a number of distributions, but is
easily extended to be able to generate random numbers from almost any of the
standard distributions.

The following table shows the classes included in the package.


CLASS STRUCTURE:

ExtReal.......................... "Extended real numbers"

Random........................... Uniform random number generator
 |
 +---Constant.................... Return a constant
 |
 +---PosGen...................... Used by PosGenX etc
 |    |
 |    +---PosGenX................ Positive random #s from decreasing density
 |    |
 |    +---Exponential............ Negative exponential rng
 |    |
 |    +---Gamma1................. Used by Gamma (shape parameter < 1)
 |    |
 |    +---SymGen................. Used by SymGenX etc
 |         |
 |         +---SymGenX........... Random #s from symmetric density
 |         |
 |         +---Cauchy............ Cauchy random number generator
 |         |
 |         +---Normal............ Standard normal random number generator
 |              |
 |              +---ChiSq1....... Used by ChiSq (one df)
 |
 +---AsymGen..................... Used by AsymGenX etc
 |    |
 |    +---AsymGenX............... Random #s from asymmetric density
 |    |
 |    +---Poisson1............... Used by Poisson (mean > 10)
 |    |
 |    +---Gamma2................. Used by Gamma (shape parameter > 1)
 |
 +---ChiSq....................... Non-central chi-squared rng
 |
 +---Gamma....................... Gamma random number generator
 |
 +---DiscreteGen................. Discrete random number generator
 |
 +---Poisson2.................... Used by Poisson (mean <= 10)
 |
 +---Poisson..................... Poisson random number generator
 |
 +---SumRandom................... Sum of random numbers
 |
 +---MixedRandom................. Mixture of random numbers
 |
 +---AddedRandom................. Used by SumRandom and MixedRandom
 |
 +---SubtedRandom................ Used by SumRandom
 |
 +---MultedRandom................ Used by SumRandom
 |
 +---ShiftedRandom............... Used by SumRandom
 |
 +---ScaledRandom................ Used by SumRandom
 |
 +---RepeatedRandom.............. Used by SumRandom
 |
 +---SelectedRandom.............. Used by MixedRandom


Each of the classes includes the following functions:

void Set(double) ................ Set seed (inherited from base class)
double Get() .................... Get seed (inherited from base class)
real Next() ..................... Get a new random number
char* Name() .................... Name of the distribution
ExtReal Mean() .................. Mean of the distribution
ExtReal Variance() .............. Variance of the distribution

plus the class constructor(s).

Set(double) must be called for just one class with an argument between 0 and 1
to set up the base random number generator before Next() is called in any
class.

The package can be run in single precision or double precision; "real" must be
defined as "float" or "double" by typedef statement. The Mean() and Variance()
functions use "extended reals" so the ExtReal package must be included.


NOW DESCRIBING THE INDIVIDUAL CLASSES AND SKETCHING THE METHODS THEY USE:

Random:
------
This is the basic uniform random number generator, used to drive all the
others. The Lewis-Goodman-Miller algorithm is used followed by Marsaglia
mixing using the C random number generator. While not perfect, and probably
superceded, the LGM generator has given me acceptable results in a wide
variety of simulations. Ideally the basic generator should be recoded in
assembly language to give the maximum speed to all the generators in this
package.

The constructor has no parameters. For example

Random U;
for (int i=0; i<100; i++) cout << U.Next() << "\n";


Constant:
--------
This returns a constant. The constructor takes one "real" parameter; the value
of the constant to be returned.


PosGen:
------
PosGen is not used directly. It is used as a base class for generating a
random number from an arbitrary probability density p(x). p(x) must be
non-zero only for x>=0, be monotonically decreasing for x>=0, and be finite.
For example, p(x) could be exp(-x) for x>=0.

The method is to cover the density in a set of rectangles of equal area as in
the diagram (indicated by ---).

|
x
|xx------
|  xx    |
|    xxx |
|.......xxx---------
|        | xxxx     |
|        |     xxxx |
|        |.........xxxxx------------
|        |          |   xxxxx       |
|        |          |        xxxxxx |
|        |          |..............xxxxxx----------------------
|        |          |               |    xxxxxxx               |
|        |          |               |           xxxxxxx        |
|        |          |               |                  xxxxxxxx|
+===========================================================================

The numbers are generated by generating a pair of numbers uniformly
distributed over these rectangles and then accepting the x coordinate as the
next random number if the pair corresponds to a point below the density
function. The acceptance can be done in two stages, the first being whether he
number is below the dotted line. This means that the density function need be
checked only very occasionally and on the average only just over 3 uniform
random numbers are required for each of the random numbers produced by this
generator.

See PosGenX or Exponential for the method of deriving a class to generate
random numbers from a given distribution.


PosGenX:
-------
This uses an arbitrary density satisfying the previous conditions to generate
random numbers from that density. Suppose

real pdf(real)

is the density. Then use pdf as the argument of the constructor. For example

PosGenX P(pdf);
for (int i=0; i<100; i++) cout << P.Next() << "\n";


Exponential:
-----------
This generates random numbers with density exp(-x) for x>=0. The constructor
takes no arguments.


Gamma1:
------
This generates random numbers from a gamma distribution with shape parameter
alpha<1. Because the density is infinite at x=0 a power transform is required
so this random number generator may not be as efficient as the others
described here. A better method would handle the immediate neighbourhood of
the origin separately and use PosGen directly for the rest. The constructor
takes real alpha as an argument.


SymGen:
------
SymGen is a modification of PosGen for unimodal distributions symmetric about
the origin, such as the standard normal.


SymGenX:
-------
This corresponds to PosGenX for symmetric distributions.


Cauchy:
------
Generates random numbers from a standard Cauchy distribution. The constructor
takes no parameters.


Normal:
------
Generates standard normal random numbers. The constructor has no arguments.
This class has been slightly modified to ensure only one copy of the arrays
generated by the constructor exist at any given time. That is, if the
constructor is called twice (before the destructor is called) only one copy of
the arrays is generated.


ChiSq1:
------
Non-central chi-squared with one degree of freedom. Used as part of ChiSq.


AsymGen:
-------
A general random number generator for unimodal distributions following the
method used by PosGen. The constructor takes one argument: the location of the
mode of the distribution.


AsymGenX:
--------
Corresponds to PosGenX. The arguments of the constructor are the name of the
density function and the location of the mode.

real pdf(real);
real mode;
.....
AsymGenX X(pdf,mode);
for (int i=0; i<100; i++) cout << X.Next() << "\n";


Poisson1:
--------
Poisson distribution; derived from AsymGen. The constructor takes the mean as
the argument. Used by Poisson for values of the mean greater than 10.


Gamma2:
------
Gamma distribution for the shape parameter, alpha, greater than 1. The
constructor takes alpha as the argument. Used by Gamma.


ChiSq:
-----
Non-Central chi-squared distribution. The method uses ChiSq1 to generate the
non-central part and Gamma2 or Exponential to generate the central part. The
constructor takes as arguments the number of degrees of freedom (>=1) and the
non-centrality parameter (omit if zero).


Gamma:
-----
Gamma distribution. The constructor takes the shape parameter as argument.
Uses Gamma1, Gamma2 or Exponential.


DiscreteGen:
-----------
This is for generating random numbers taking just a finite number of values.
There are two alternative forms of the constructor:

DiscreteGen D(n,prob);
DiscreteGen D(n,prob,val);

where n is an integer giving the number of values, prob a real array of length
n giving the probabilities and val a real array of length n giving the set of
values that are generated. If val is omitted the values are 0,1,...,n-1.

The method requires two uniform random numbers for each number it produces.
This method is described by Kronmal and Peterson, American Statistician, 1979,
Vol 33, No 4, pp214-218.


Poisson2:
--------
Poisson distribution with mean less than or equal to 10. Uses DiscreteGen.
Constructor takes the mean as its argument.


Poisson:
-------
Poisson distribution: uses Poisson1 or Poisson2. Constructor takes the mean as
its argument.


SumRandom:
---------
This is for building a random number generator as a linear or multplicative
combination of existing random number generators. Suppose rv1, rv2, rv3, rv4
are random number generators defined with constructors given above and r1, r2,
r0 are reals and i1, i3 are integers.

Then the generator S defined by

SumRandom S = rv1(i1)*r1 - rv2*r2 + rv3(i3)*rv4 + r0;

has the obvious meaning. rv1(i1) means that the sum of i1 independent values
from of rv1 should be used. Note that rv1*rv1 means the product of two
independent numbers generated from rv1.

Remember that this becomes a very inefficient method if the number of terms or
copies is large.


MixedRandom:
-----------
This is for mixtures of distributions. Suppose rv1, rv2, rv3 are random number
generators and p1, p2, p3 are reals summing to 1.

Then the generator M defined by

MixedRandom M = rv1(p1) + rv2(p2) + rv3(p3);

produces a random number generator with selects its next random number from
rv1 with probability p1, rv2 with probability p2, rv3 with probability p3.

Alternatively one can use the constructor

MixedRandom M(n, prob, rv);

where n is the number of distributions in the mixture, prob the real array of
probabilities, rv an array of pointers to random variables.


AddedRandom, SubtedRandom, MultedRandom, ShiftedRandom, ScaledRandom,
-----------  ------------  ------------  -------------  ------------
RepeatedRandom, SelectedRandom:
--------------  --------------
These are used by SumRandom and MixedRandom.


GENERATING RANDOM NUMBERS FROM OTHER DISTRIBUTIONS:

Continuous finite unimodal density (no parameters): I assume you have a
program to calculate the density. Use PosGenX, SymGenX or AsymGen.
Alternatively derive a new class from PosGen, SymGen or AsymGen. Cauchy and
Exponential are examples.

Continuous finite unimodal density (with parameters): Derive a new class from
PosGen, SymGen or AsymGen. Gamma2 is an example.

Transformation of existing random number generator: Derive a new class from
the existing class. ChiSq1 is an example.

Transformation of several random numbers: Define a new class derived from
Random and generate the new random number from the existing random number
generators. ChiSq is a rather complex example.

Density with infinite singularity: This may sometimes be handled by
transforming a random variable generated by PosGen, SymGen or AsymGen. Gamma1
is an example.

Distribution with several modes: Breakdown into a mixture of unimodal
distributions.

Linear or quadratic combination of existing random numbers: Use SumRandom.

Mixture of existing random numbers: Use MixedRandom.

Discrete distribution (< 100 possible values): Use DiscreteGen. DiscreteGen
and Poisson2 are examples.

Discrete distribution (many possible values): Use PosGen, SymGen or AsymGen.
Poisson1 is an example.


MAKEFILE TO COMPILE UNDER ZORTECH:

Assume that main() is in file tryrand.cxx.


C = ztc -f -c -dZORTECH $*.cxx
D = ztc -f -c -dZORTECH $*.cxx -o

OBJ = tryrand.obj random.obj extreal.obj

makefile:     tryrand.exe

tryrand.exe:  $(OBJ)
              ztc -f $**

random.obj:   random.cxx extreal.hxx random.hxx
              $D

extreal.obj:  extreal.cxx extreal.hxx
              $D

tryrand.obj:  tryrand.cxx
              $D


FILES INCLUDED IN THIS PACKAGE:

random.tex         this file
random.hxx         include file for random.cxx
random.cxx         main code file
extreal.hxx        include file for extreal.cxx
extreal.cxx        code file for extreal package
boolean.hxx        definition of boolean variable


// random.hxx ------------------------------------------------------------


#ifndef RANDOM_LIB
#define RANDOM_LIB 0

#ifndef ZORTECH                  // for systems that won't accept long names
#define AddedRandom AdRan
#define MultedRandom MuRan
#define SubtedRandom SuRan
#define ScaledRandom ScRan
#define ShiftedRandom ShRan
#define RepeatedRandom ReRan
#define SelectedRandom SeRan
#define SumRandom SumRan
#define MixedRandom MiRan
#endif



/******************** utilities and definitions *************************/

#include <boolean.hxx>
#include <extreal.hxx>

typedef real (*PDF)(real);                // probability density function

/***************** uniform random number generator **********************/

class AddedRandom;                        // defined later
class MultedRandom;
class SubtedRandom;
class ScaledRandom;
class ShiftedRandom;
class RepeatedRandom;
class SelectedRandom;

class Random                              // uniform random number generator
{
   static double seed;                    // seed
   static real Buffer[128];               // for mixing random numbers
   real Raw();                            // unmixed random numbers

protected:
   virtual void Error(char*);             // error handler
   virtual void ErrorNoSpace();           // nospace message

public:
   void Set(double s);                    // set seed (0 < seed < 1)
   double Get();                          // get seed
   virtual real Next();                   // get new value
   virtual char* Name();                  // identification
   virtual real Density(real);            // used by PosGen & Asymgen
   Random() {}                            // do nothing
   virtual ~Random() {}                   // make destructors virtual
   AddedRandom operator+(Random&);        // for making combinations
   MultedRandom operator*(Random&);
   SubtedRandom operator-(Random&);
   ShiftedRandom operator+(real);
   ShiftedRandom operator-(real);
   ScaledRandom operator*(real);
   RepeatedRandom operator()(int);
   SelectedRandom operator()(double);     // used by MixedRandom
   void operator=(Random);                // = (not allowed)
   virtual ExtReal Mean() { return 0.5; } // mean of distribution
   virtual ExtReal Variance() { return 1.0/12.0; }
                                          // variance of distribution
   virtual int nelems() { return 1; }     // used by MixedRandom
   virtual void load(int*,real*,Random**);
};

/************************** return constant ******************************/

class Constant : public Random
{
   real value;                            // value to be returned

public:
   char* Name();                          // identification
   Constant(real v) { value=v; }          // set value
   real Next() { return value; }
   ExtReal Mean() { return value; }
   ExtReal Variance() { return 0.0; }
};

/***************** positive random number generator **********************/

class PosGen : public Random              // generate positive rv
{
protected:
   real xi, *sx, *sfx;
   BOOL NotReady;
   void Build(BOOL);                      // called on first call to Next

public:
   char* Name();                          // identification
   PosGen();                              // constructor
   ~PosGen();                             // destructor
   real Next();                           // to get a single new value
   ExtReal Mean() { return Missing; }
   ExtReal Variance() { return Missing; }
};

/***************** symmetric random number generator **********************/

class SymGen : public PosGen              // generate symmetric rv
{

public:
   char* Name();                          // identification
   real Next();                           // to get a single new value
};

/***************** normal random number generator **********************/

class Normal : public SymGen              // generate standard normal rv
{
   static real Nxi, *Nsx, *Nsfx;          // so we need initialise only once
   static long count;                     // assume initialised to 0
public:
   char* Name();                          // identification
   Normal();
   ~Normal();
   real Density(real);                    // normal density function
   ExtReal Mean() { return 0.0; }
   ExtReal Variance() { return 1.0; }
};

/************ chi-square random number generator **********************/

class ChiSq : public Random               // generate non-central chi-sq rv
{
   Random* c1;                            // pointers to generators
   Random* c2;                            // pointers to generators
   int version;                           // indicates method of generation
   real mean, var;

public:
   char* Name();                          // identification
   ChiSq(int, real=0.0);                  // df and non-centrality parameter
   ~ChiSq();
   ExtReal Mean() { return mean; }
   ExtReal Variance() { return var; }
   real Next();
};

/***************** Cauchy random number generator **********************/

class Cauchy : public SymGen              // generate standard cauchy rv
{
public:
   char* Name();                          // identification
   real Density(real);                    // Cauchy density function
   ExtReal Mean() { return Indefinite; }
   ExtReal Variance() { return PlusInfinity; }
};

/***************** Exponential random number generator **********************/

class Exponential : public PosGen         // generate standard cauchy rv
{
public:
   char* Name();                          // identification
   real Density(real);                    // Cauchy density function
   ExtReal Mean() { return 1.0; }
   ExtReal Variance() { return 1.0; }
};

/***************** asymmetric random number generator **********************/

class AsymGen : public Random             // generate asymmetric rv
{
   real xi, *sx, *sfx, mode; int ic;
   BOOL NotReady;
   void Build();                          // called on first call to Next

public:
   char* Name();                          // identification
   AsymGen(real);                         // constructor (real=mode)
   ~AsymGen();                            // destructor
   real Next();                           // to get a single new value
   ExtReal Mean() { return Missing; }
   ExtReal Variance() { return Missing; }
};

/***************** Gamma random number generator **********************/

class Gamma : public Random               // generate gamma rv
{
   Random* method;
public:
   char* Name();                          // identification
   Gamma(real);                           // constructor (real=shape)
   ~Gamma() { delete method; }
   real Next() { return method->Next(); }
   ExtReal Mean() { return method->Mean(); }
   ExtReal Variance() { return method->Variance(); }
};

/***************** Generators with pointers to pdf ***********************/

class PosGenX : public PosGen
{
   PDF f;
public:
   char* Name();                          // identification
   PosGenX(PDF fx) { f=fx; }
   real Density(real x) { return (*f)(x); }
};

class SymGenX : public SymGen
{
   PDF f;
public:
   char* Name();                          // identification
   SymGenX(PDF fx) { f=fx; }
   real Density(real x) { return (*f)(x); }
};

class AsymGenX : public AsymGen
{
   PDF f;
public:
   char* Name();                          // identification
   AsymGenX(PDF fx, real mode) : (mode) { f=fx; }
   real Density(real x) { return (*f)(x); }
};

/***************** discrete random number generator **********************/

class DiscreteGen : public Random         // discrete random generator
{
   real *p; int *ialt; int n; real *val;
   void Gen(int, real*);                  // called by constructors
   real mean, var;                        // calculated by the constructor

public:
   char* Name();                          // identification
   DiscreteGen(int,real*);                // constructor
   DiscreteGen(int,real*,real*);          // constructor
   ~DiscreteGen();                        // destructor
   real Next();                           // new single value
   ExtReal Mean() { return mean; }
   ExtReal Variance() { return var; }
};

/****************** Poisson random number generator *******************/

class Poisson : public Random            // generate poisson rv
{
   Random* method;
public:
   char* Name();                          // identification
   Poisson(real);                         // constructor (real=mean)
   ~Poisson() { delete method; }
   real Next() { return method->Next(); }
   ExtReal Mean() { return method->Mean(); }
   ExtReal Variance() { return method->Variance(); }
};

/************************* sum of random numbers ************************/

class ScaledRandom : public Random
{
   friend Random;
   Random* rv; real s;
   ScaledRandom(ScaledRandom& r) { rv=r.rv; s=r.s; }
   ScaledRandom() {}
public:
   real Next();
   ExtReal Mean();
   ExtReal Variance();
};

class ShiftedRandom : public Random
{
   friend Random;
   Random* rv; real s;
   ShiftedRandom(ShiftedRandom& r) { rv=r.rv; s=r.s; }
   ShiftedRandom() {}
public:
   real Next();
   ExtReal Mean();
   ExtReal Variance();
};

class RepeatedRandom : public Random
{
   friend Random;
   Random* rv; int n;
   RepeatedRandom(RepeatedRandom& r) { rv=r.rv; n=r.n; }
   RepeatedRandom() {}
public:
   real Next();
   ExtReal Mean();
   ExtReal Variance();
};

class AddedRandom : public Random
{
   friend Random;
   Random* rv1; Random* rv2;
   AddedRandom(AddedRandom& r) { rv1=r.rv1; rv2=r.rv2; }
   AddedRandom() {}
public:
   real Next();
   ExtReal Mean();
   ExtReal Variance();
   int nelems();
   void load(int*,real*,Random**);
};

class MultedRandom : public Random
{
   friend Random;
   Random* rv1; Random* rv2;
   MultedRandom(MultedRandom& r) { rv1=r.rv1; rv2=r.rv2; }
   MultedRandom() {}
public:
   real Next();
   ExtReal Mean();
   ExtReal Variance();
};

class SubtedRandom : public Random
{
   friend Random;
   Random* rv1; Random* rv2;
   SubtedRandom(SubtedRandom& r) { rv1=r.rv1; rv2=r.rv2; }
   SubtedRandom() {}
public:
   real Next();
   ExtReal Mean();
   ExtReal Variance();
};

class SumRandom : public Random           // sum of random variables
                                          // version user can use
{
   Random* rv;
public:
   char* Name();                          // identification
   SumRandom(AddedRandom& r) { rv=&r; }
   SumRandom(MultedRandom& r) { rv=&r; }
   SumRandom(SubtedRandom& r) { rv=&r; }
   SumRandom(ShiftedRandom& r) { rv=&r; }
   SumRandom(ScaledRandom& r) { rv=&r; }
   SumRandom(RepeatedRandom& r) { rv=&r; }
   real Next() { return rv->Next(); }
   ExtReal Mean() { return rv->Mean(); }
   ExtReal Variance() { return rv->Variance(); }
};

/********************* mixtures of random numbers ************************/

class SelectedRandom : public Random
{
   friend Random;
   Random* rv; real p;
   SelectedRandom(SelectedRandom& r) { rv=r.rv; p=r.p; }
   SelectedRandom() {}
public:
   void load(int*,real*,Random**);
   real Next();
};

class MixedRandom : public Random         // mixtures of random numbers
{
   int n;                                 // number of components
   DiscreteGen* dg;                       // discrete mixing distribution
   Random** rv;                           // array of pointers to rvs
   ExtReal mean, var;
   void Build(real*);                     // used by constructors

public:
   char* Name();                          // identification
   MixedRandom(int, real*, Random**);
   MixedRandom(AddedRandom&);
   ~MixedRandom();
   real Next();
   ExtReal Mean() { return mean; }
   ExtReal Variance() { return var; }
};

#endif

// random.cxx -----------------------------------------------------------


typedef double real;                         // all floating point double

//#define MONITOR = 0                // monitor constructors and destructors

#ifdef ZORTECH
   #include <math.h>
   #include <stdlib.h>
   #include <stream.hpp>
#else
   #include <math.hxx>
   #include <xstream.hxx>
   int rand();
#endif

#include <random.hxx>

/************ chi-square-1 random number generator **********************/

class ChiSq1 : public Normal              // generate non-central chi-square
                                          // rv with 1 degree of freedom
{
   real deltasq;                          // non-centrality parameter
   real delta;

public:
   ChiSq1(real);                          // non-centrality parameter
   ExtReal Mean() { return 1.0+deltasq; }
   ExtReal Variance() { return 2.0+4.0*deltasq; }
   real Next();
};

/************ Poisson random number generator, larger mu ****************/

class Poisson1 : public AsymGen           // generate poisson rv, large mu
{
   real mu, ln_mu;
public:
   Poisson1(real);                        // constructor (real=mean)
   real Density(real);                    // Poisson density function
   real Next();                           // truncate value from AsymGen
   ExtReal Mean() { return mu; }
   ExtReal Variance() { return mu; }
};

/*********** Gamma random number generator, alpha <= 1 *****************/

class Gamma1 : public PosGen              // generate gamma rv
                                          // shape parameter <= 1
{
   real ln_gam, ralpha, alpha;
public:
   Gamma1(real);                          // constructor (real=shape)
   real Density(real);                    // gamma density function
   real Next();                           // carries out power transform
   ExtReal Mean() { return alpha; }
   ExtReal Variance() { return alpha; }
};

/*********** Gamma random number generator, alpha > 1 ******************/

class Gamma2 : public AsymGen             // generate gamma rv
                                          // shape parameter > 1
{
   real alpha, ln_gam;
public:
   Gamma2(real);                          // constructor (real=shape)
   real Density(real);                    // gamma density function
   ExtReal Mean() { return alpha; }
   ExtReal Variance() { return alpha; }
};

/************ Poisson random number generator, mu <= 10 ****************/

class Poisson2 : public Random            // generate poisson rv, large mu
{
   DiscreteGen* dg;
public:
   Poisson2(real);                        // constructor (real=mean)
   ~Poisson2();
   real Next() { return dg->Next(); }
   ExtReal Mean() { return dg->Mean(); }
   ExtReal Variance() { return dg->Variance(); }
};

/***************************** utilities ******************************/

real ln_gamma(real);                      // log gamma function

overload Square;
inline real Square(real x) { return x*x; }
inline ExtReal Square(ExtReal& x) { return x*x; }

/************************** end of definitions ************************/

real Random::Raw()                           // get new uniform random number
{
   // m = 2147483647 = 2^31 - 1; a = 16807;
   // 127773 = m div a; 2836 = m mod a
   long iseed = (long)seed;
   long hi = iseed / 127773;                 // integer division
   long lo = iseed - hi * 127773;            // modulo
   iseed = 16807 * lo - 2836 * hi;
   if (iseed <= 0) iseed += 2147483647;
   seed = (double)iseed; return seed*4.656612875e-10;
}

real Random::Density(real)
{ Error("density function not defined"); return 0.0; }

real Random::Next()                          // get new mixed random number
{
   int i = rand()/256;                       // 0 <= i < 128
   real f = Buffer[i]; Buffer[i] = Raw();
   return f;
}

double Random::Get()                         // get random number seed
{ return seed*4.656612875e-10; }

void Random::Set(double s)                   // set random number seed
                                             // s must be between 0 and 1
{
   if (s>=1.0 || s<=0.0) Error("seed out of range");
   seed = (long)(s/4.656612875e-10);
   for (int i = 0; i<128; i++) Buffer[i] = Raw();
}

void Random::operator=(Random) { Error("'=' not allowed"); }

void Random::Error(char* message)
{
   cerr << "\nError in " << this->Name() << ": " << message << "\n";
   exit(1);
}

void Random::ErrorNoSpace()
{
   cerr << "\nError in " << this->Name() << ": no space\n";
   exit(1);
}


PosGen::PosGen()                             // Constructor
{
   #ifdef MONITOR
      cout << "constructing PosGen\n";
   #endif
   NotReady=TRUE;
}

PosGen::~PosGen()
{
   if (!NotReady)
   {
      #ifdef MONITOR
         cout << "freeing PosGen arrays\n";
      #endif
      delete [60] sx; delete [60] sfx;
   }
   #ifdef MONITOR
      cout << "destructing PosGen\n";
   #endif
}

void PosGen::Build(BOOL sym)                 // set up arrays
{
   #ifdef MONITOR
      cout << "building PosGen arrays\n";
   #endif
   NotReady=FALSE;
   sx=new real[60]; sfx=new real[60];
   if (!sx || !sfx) ErrorNoSpace();
   real sxi=0.0; real inc = sym ? 0.01 : 0.02;
   for (int i=0; i<60; i++)
   {
      sx[i]=sxi; real f1=Density(sxi); sfx[i]=f1;
      if (f1<=0.0) goto L20;
      sxi+=inc/f1;
   }
   Error("area too large");
L20:
   if (i<50) Error("area too small");
   xi = sym ? 2*i : i;
   return;
}

real PosGen::Next()
{
   real ak,y; int ir;
   if (NotReady) Build(FALSE);
   do {
      real r1=Random::Next();
      ir = (int)(r1*xi); real sxi=sx[ir];
      ak=sxi+(sx[ir+1]-sxi)*Random::Next();
      y=sfx[ir]*Random::Next();
   } while ( y>=sfx[ir+1] && y>=Density(ak) );
   return ak;
}

real SymGen::Next()
{
   real s,ak,y; int ir;
   if (NotReady) Build(TRUE);
   do {
      s=1.0;
      real r1=Random::Next();
      if (r1>0.5) { s=-1.0; r1=1.0-r1; }
      ir = (int)(r1*xi); real sxi=sx[ir];
      ak=sxi+(sx[ir+1]-sxi)*Random::Next();
      y=sfx[ir]*Random::Next();
   } while ( y>=sfx[ir+1] && y>=Density(ak) );
   return s*ak;
}

AsymGen::AsymGen(real modex)                 // Constructor
{
   #ifdef MONITOR
      cout << "constructing AsymGen\n";
   #endif
   mode=modex; NotReady=TRUE;
}

void AsymGen::Build()                        // set up arrays
{
   #ifdef MONITOR
      cout << "building AsymGen arrays\n";
   #endif
   NotReady=FALSE;
   sx=new real[121]; sfx=new real[121];
   if (!sx || !sfx)  ErrorNoSpace();
   real sxi=mode;
   for (int i=0; i<120; i++)
   {
      sx[i]=sxi; real f1=Density(sxi); sfx[i]=f1;
      if (f1<=0.0) goto L20;
      sxi+=0.01/f1;
   }
   Error("area too large (a)");
L20:
   ic=i-1; sx[120]=sxi; sfx[120]=0.0;
   sxi=mode;
   for (; i<120; i++)
   {
      sx[i]=sxi; real f1=Density(sxi); sfx[i]=f1;
      if (f1<=0.0) goto L30;
      sxi-=0.01/f1;
   }
   Error("area too large (b)");
L30:
   if (i<100)  Error("area too small");
   xi=i;
   return;
}

real AsymGen::Next()
{
   real ak,y; int ir1;
   if (NotReady) Build();
   do {
      real r1=Random::Next();
      int ir=(int)(r1*xi); real sxi=sx[ir];
      ir1 = (ir==ic) ? 120 : ir+1;
      ak=sxi+(sx[ir1]-sxi)*Random::Next();
      y=sfx[ir]*Random::Next();
   } while ( y>=sfx[ir1] && y>=Density(ak) );
   return ak;
}

AsymGen::~AsymGen()
{
   if (!NotReady)
   {
      #ifdef MONITOR
         cout << "freeing AsymGen arrays\n";
      #endif
      delete [121] sx; delete [121] sfx;
   }
   #ifdef MONITOR
      cout << "destructing AsymGen\n";
   #endif
}

Normal::Normal()
{
   if (count) { NotReady=FALSE; xi=Nxi; sx=Nsx; sfx=Nsfx; }
   else { Build(TRUE); Nxi=xi; Nsx=sx; Nsfx=sfx; }
   count++;
}

Normal::~Normal()
{
   count--;
   if (count) NotReady=TRUE;                     // disable freeing arrays
}

real Normal::Density(real x)                     // normal density
{ return (fabs(x)>8.0) ? 0 : 0.398942280 * exp(-x*x / 2); }

ChiSq1::ChiSq1(real d) : ()                      // chisquare with 1 df
{ deltasq=d; delta=sqrt(d); }

real ChiSq1::Next()
{ real z=Normal::Next()+delta; return z*z; }

ChiSq::ChiSq(int df, real noncen)
{
  if (df<=0 || noncen<0.0) Error("illegal parameters");
  else if (df==1) { version=0; c1=new ChiSq1(noncen); }
  else if (noncen==0.0)
  {
     if (df==2) { version=1; c1=new Exponential(); }
     else { version=2; c1=new Gamma2(0.5*df); }
  }
  else if (df==2) { version=3; c1=new ChiSq1(noncen/2.0); }
  else if (df==3) { version=4; c1=new Exponential(); c2=new ChiSq1(noncen); }
  else { version=5; c1=new Gamma2(0.5*(df-1)); c2=new ChiSq1(noncen); }
  if (!c1 || (version>3 && !c2)) ErrorNoSpace();
  mean=df+noncen; var=2*df+4.0*noncen;
}

ChiSq::~ChiSq() { delete c1; if (version>3) delete c2; }

real ChiSq::Next()
{
   switch(version)
   {
   case 0: return c1->Next();
   case 1: case 2: return c1->Next()*2.0;
   case 3: return c1->Next() + c1->Next();
   case 4: case 5: return c1->Next()*2.0 + c2->Next();
   }
}

real Cauchy::Density(real x)                     // Cauchy density function
{ return (fabs(x)>1.0e15) ? 0 : 0.31830988618 / (1.0+x*x); }

Poisson1::Poisson1(real mux) : (mux)             // Constructor
{ mu=mux; ln_mu=log(mu); }

real Poisson1::Density(real x)                   // Poisson density function
{
   if (x<0.0) return 0.0;
   double ix = floor(x);                         // use integer part
   double l = ln_mu * ix - mu - ln_gamma(1.0+ix);
   return  (l < -30.0) ? 0.0 : exp(l);
}

real Poisson1::Next() { return floor(AsymGen::Next()); }

Poisson2::Poisson2(real mux)
{
   real probs[40];
   probs[0]=exp(-mux);
   for (int i=1; i<40; i++) probs[i]=probs[i-1]*mux/i;
   dg=new DiscreteGen(40,probs);
   if (!dg) ErrorNoSpace();
}

Poisson2::~Poisson2() { delete dg; }

real Exponential::Density(real x)                // Negative exponential
{ return  (x > 30.0) ? 0.0 : exp(-x); }

Poisson::Poisson(real mu)
{
   if (mu<=10.0) method = new Poisson2(mu);
   else method = new Poisson1(mu);
   if (!method) ErrorNoSpace();
}

Gamma1::Gamma1(real alphax)                      // constructor (real=shape)
{ ralpha=1.0/alphax; ln_gam=ln_gamma(alphax+1.0); alpha=alphax; }

real Gamma1::Density(real x)                     // density function for
{                                                // transformed gamma
   double l = - pow(x,ralpha) - ln_gam;
   return  (l < -30.0) ? 0.0 : exp(l);
}

real Gamma1::Next()                               // transform variable
{ return pow(PosGen::Next(),ralpha); }

Gamma2::Gamma2(real alphax) : (alphax-1.0)       // constructor (real=shape)
{ alpha=alphax; ln_gam=ln_gamma(alpha); }

real Gamma2::Density(real x)                     // gamma density function
{
   if (x<=0.0) return 0.0;
   double l = (alpha-1.0)*log(x) - x - ln_gam;
   return  (l < -30.0) ? 0.0 : exp(l);
}

Gamma::Gamma(real alpha)                         // general gamma generator
{
   if (alpha<1.0) method = new Gamma1(alpha);
   else if (alpha==1.0) method = new Exponential();
   else method = new Gamma2(alpha);
   if (!method)  ErrorNoSpace();
}

DiscreteGen::DiscreteGen(int n1, real* prob)     // discrete generator
                                                 // values on 0,...,n1-1
{
   #ifdef MONITOR
      cout << "constructing DiscreteGen\n";
   #endif
   Gen(n1, prob); val=0;
   mean=0.0; var=0.0;
   { for (int i=0; i<n; i++) mean = mean + i*prob[i]; }
   { for (int i=0; i<n; i++) var = var + Square(i-mean) * prob[i]; }
}

DiscreteGen::DiscreteGen(int n1, real* prob, real* val1)
                                                 // discrete generator
                                                 // values on *val
{
   #ifdef MONITOR
      cout << "constructing DiscreteGen\n";
   #endif
   Gen(n1, prob); val = new real[n1];
   if (!val)  ErrorNoSpace();
   for (int i=0; i<n1; i++) val[i]=val1[i];
   mean=0.0; var=0.0;
   { for (int i=0; i<n; i++) mean = mean + val[i]*prob[i]; }
   { for (int i=0; i<n; i++) var = var + Square(val[i]-mean)*prob[i]; }
}


void DiscreteGen::Gen(int n1, real* prob)
{
   n=n1;                                         // number of values
   p=new real[n]; ialt=new int[n];
   if (!p || !ialt)  ErrorNoSpace();
   real rn = 1.0/n; real px; int i;
   for (i=0; i<n; i++) { p[i]=0.0; ialt[i]=-1; }
   for (i=0; i<n; i++)
   {
      real pmin=1.0; real pmax=-1.0; int jmin=-1; int jmax=-1;
      for (int j=0; j<n; j++)
      {
         if (ialt[j]<0)
         {
            px=prob[j]-p[j];
            if (pmax<=px) { pmax=px; jmax=j; }
            if (pmin>=px) { pmin=px; jmin=j; }
         }
      }
      if ((jmax<0) || (jmin<0)) Error("method fails");
      ialt[jmin]=jmax; px=rn-pmin; p[jmax]+=px; px*=n; p[jmin]=px;
      if ((px>1.00001)||(px<-.00001)) Error("probs don't add to 1 (a)");
   }
   if (px>0.00001) Error("probs don't add to 1 (b)");
}

DiscreteGen::~DiscreteGen()
{
   delete [n] p; delete [n] ialt; if (!val) delete [n] val;
   #ifdef MONITOR
      cout << "destructing DiscreteGen\n";
   #endif
}

real DiscreteGen::Next()                  // Next discrete random variable
{
   int i = (int)(n*Random::Next()); if (Random::Next()<p[i]) i=ialt[i];
   return val ? val[i] : (real)i;
}

real ln_gamma(real xx)
{
   // log gamma function from numerical recipies in C

   if (xx<1.0)                           // Use reflection formula
   {
      const double pi=3.14159265359;
      double piz = pi * (1.0-xx);
      return log(piz/sin(piz))-ln_gamma(2.0-xx);
   }
   else
   {
      static double cof[6]={76.18009173,-86.50532033,24.01409822,
         -1.231739516,0.120858003e-2,-0.536382e-5};

      double x=xx-1.0; double tmp=x+5.5;
      tmp -= (x+0.5)*log(tmp); double ser=1.0;
      for (int j=0;j<=5;j++) { x += 1.0; ser += cof[j]/x; }
      return -tmp+log(2.50662827465*ser);
   }
}

real ScaledRandom::Next() { return rv->Next()*s; }

ExtReal ScaledRandom::Mean() { return rv->Mean()*s; }

ExtReal ScaledRandom::Variance() { return rv->Variance()*(s*s); }

real ShiftedRandom::Next() { return rv->Next()+s; }

ExtReal ShiftedRandom::Mean() { return rv->Mean()+s; }

ExtReal ShiftedRandom::Variance() { return rv->Variance(); }

ExtReal RepeatedRandom::Mean() { return rv->Mean()*(real)n; }

ExtReal RepeatedRandom::Variance() { return rv->Variance()*(real)n; }

RepeatedRandom Random::operator()(int n)
{ RepeatedRandom r; r.rv=this; r.n=n; return r; }

ShiftedRandom Random::operator+(real s)
{ ShiftedRandom r; r.rv=this; r.s=s; return r; }

ShiftedRandom Random::operator-(real s)
{ ShiftedRandom r; r.rv=this; r.s=-s; return r; }

ScaledRandom Random::operator*(real s)
{ ScaledRandom r; r.rv=this; r.s=s; return r; }

AddedRandom Random::operator+(Random& rv2)
{ AddedRandom r; r.rv1=this; r.rv2=&rv2; return r; }

MultedRandom Random::operator*(Random& rv2)
{ MultedRandom r; r.rv1=this; r.rv2=&rv2; return r; }

SubtedRandom Random::operator-(Random& rv2)
{ SubtedRandom r; r.rv1=this; r.rv2=&rv2; return r; }

real AddedRandom::Next() { return rv1->Next() + rv2->Next() ; }

ExtReal AddedRandom::Mean() { return rv1->Mean() + rv2->Mean() ; }

ExtReal AddedRandom::Variance() { return rv1->Variance() + rv2->Variance() ; }

int AddedRandom::nelems() { return rv1->nelems() + rv2->nelems(); }

void AddedRandom::load(int* i, real* probs, Random** rvx)
{ rv1->load(i, probs, rvx); rv2->load(i, probs, rvx); }

real SubtedRandom::Next() { return rv1->Next() - rv2->Next() ; }

ExtReal SubtedRandom::Mean() { return rv1->Mean() - rv2->Mean() ; }

ExtReal SubtedRandom::Variance() { return rv1->Variance() + rv2->Variance() ; }

real MultedRandom::Next() { return rv1->Next() * rv2->Next() ; }

ExtReal MultedRandom::Mean() { return rv1->Mean() * rv2->Mean() ; }

ExtReal MultedRandom::Variance()
{
   ExtReal e1 = Square(rv1->Mean()); ExtReal e2 = Square(rv2->Mean());
   ExtReal v1 = rv1->Variance(); ExtReal v2 = rv2->Variance();
   ExtReal r=v1*v2; r=r+v1*e2; r=r+e1*v2;  // sum doesn't work under Zortech
   return r;
}

void Random::load(int*,real*,Random**) { Error("Illegal combination"); }

void SelectedRandom::load(int* i, real* probs, Random** rvx)
{ probs[*i]=p; rvx[*i]=rv; (*i)++; }

real SelectedRandom::Next() { Error("Next not defined"); return 0.0; }

real RepeatedRandom::Next()
{ real sum=0.0; for (int i=0; i<n; i++) sum+=rv->Next(); return sum; }

MixedRandom::MixedRandom(int nx, real* probs, Random** rvx)
{
   n=nx;
   rv=new Random*[n]; if (!rv) ErrorNoSpace();
   for (int i=0; i<n; i++) rv[i]=rvx[i];
   Build(probs);
}

MixedRandom::MixedRandom(AddedRandom& sr)
{
   n = sr.nelems();                       // number of terms;
   real* probs = new real[n]; rv = new Random*[n];
   if (!probs || !rv) ErrorNoSpace();
   int i=0;
   sr.load(&i,probs,rv);
   Build(probs);
   delete [n] probs;
}

void MixedRandom::Build(real* probs)
{
   dg=new DiscreteGen(n,probs);
   if (!dg) ErrorNoSpace();
   mean=0.0; var=0.0;
   for (int i=0; i<n; i++) mean = mean + (rv[i])->Mean()*probs[i];
   for (i=0; i<n; i++)
   {
      ExtReal sigsq=(rv[i])->Variance();
      ExtReal mudif=(rv[i])->Mean()-mean;
      var = var + (sigsq+Square(mudif))*probs[i];
   }
}

MixedRandom::~MixedRandom() { delete [n] rv; delete dg; }

real MixedRandom::Next()
{
   int i = (int)(0.5+dg->Next());          // rounding
   return (rv[i])->Next();
}

SelectedRandom Random::operator()(double p)
{ SelectedRandom r; r.rv=this; r.p=p; return r; }


// Identification routines for each class - may be not work on all compilers

char* Random::Name()       { return "Random";      }
char* Constant::Name()     { return "Constant";    }
char* PosGen::Name()       { return "PosGen";      }
char* SymGen::Name()       { return "SymGen";      }
char* AsymGen::Name()      { return "AsymGen";     }
char* PosGenX::Name()      { return "PosGenX";     }
char* SymGenX::Name()      { return "SymGenX";     }
char* AsymGenX::Name()     { return "AsymGenX";    }
char* Normal::Name()       { return "Normal";      }
char* ChiSq::Name()        { return "ChiSq";       }
char* Cauchy::Name()       { return "Cauchy";      }
char* Exponential::Name()  { return "Exponential"; }
char* Poisson::Name()      { return "Poisson";     }
char* Gamma::Name()        { return "Gamma";       }
char* DiscreteGen::Name()  { return "DiscreteGen"; }
char* SumRandom::Name()    { return "SumRandom";   }
char* MixedRandom::Name()  { return "MixedRandom"; }

// extreal.hxx ----------------------------------------------------------


#ifndef EXTREAL_LIB
#define EXTREAL_LIB 0

/************************ extended real class ***************************/

enum EXT_REAL_CODE
   { Finite, PlusInfinity, MinusInfinity, Indefinite, Missing };

class ExtReal
{
   real value;
   EXT_REAL_CODE c;
public:
   ExtReal operator+(ExtReal&);
   ExtReal operator-(ExtReal&);
   ExtReal operator*(ExtReal&);
   ExtReal operator-();
   friend ostream& operator<<( ostream&, ExtReal& );
   ExtReal(real v) { c=Finite; value=v; }
   ExtReal(EXT_REAL_CODE cx) { c=cx; }
   ExtReal() { c = Missing; }
};

#endif

// extreal.cxx ----------------------------------------------------------


typedef double real;

#ifdef ZORTECH
   #include <stream.hpp>
#else
   #include <xstream.hxx>
#endif

#include <extreal.hxx>

ExtReal ExtReal::operator+(ExtReal& er)
{
   if (c==Finite && er.c==Finite) return ExtReal(value+er.value);
   if (c==Missing || er.c==Missing) return ExtReal(Missing);
   if (c==Indefinite || er.c==Indefinite) return ExtReal(Indefinite);
   if (c==PlusInfinity)
   {
      if (er.c==MinusInfinity) return ExtReal(Indefinite);
      return *this;
   }
   if (c==MinusInfinity)
   {
      if (er.c==PlusInfinity) return ExtReal(Indefinite);
      return *this;
   }
   return er;
}

ExtReal ExtReal::operator-(ExtReal& er)
{
   if (c==Finite && er.c==Finite) return ExtReal(value-er.value);
   if (c==Missing || er.c==Missing) return ExtReal(Missing);
   if (c==Indefinite || er.c==Indefinite) return ExtReal(Indefinite);
   if (c==PlusInfinity)
   {
      if (er.c==PlusInfinity) return ExtReal(Indefinite);
      return *this;
   }
   if (c==MinusInfinity)
   {
      if (er.c==MinusInfinity) return ExtReal(Indefinite);
      return *this;
   }
   return er;
}

ExtReal ExtReal::operator*(ExtReal& er)
{
   if (c==Finite && er.c==Finite) return ExtReal(value*er.value);
   if (c==Missing || er.c==Missing) return ExtReal(Missing);
   if (c==Indefinite || er.c==Indefinite) return ExtReal(Indefinite);
   if (c==Finite)
   {
      if (value==0.0) return ExtReal(Indefinite);
      if (value>0.0) return er;
      return (-er);
   }
   if (er.c==Finite)
   {
      if (er.value==0.0) return ExtReal(Indefinite);
      if (er.value>0.0) return *this;
      return -(*this);
   }
   if (c==PlusInfinity) return er;
   return (-er);
}

ExtReal ExtReal::operator-()
{
   switch (c)
   {
      case Finite:        return ExtReal(-value);
      case PlusInfinity:  return ExtReal(MinusInfinity);
      case MinusInfinity: return ExtReal(PlusInfinity);
      case Indefinite:    return ExtReal(Indefinite);
      case Missing:       return ExtReal(Missing);
   }
}

ostream& operator<<(ostream& os, ExtReal& er)
{
   switch (er.c)
   {
      case Finite:        os << er.value;         break;
      case PlusInfinity:  os << "plus-infinity";  break;
      case MinusInfinity: os << "minus-infinity"; break;
      case Indefinite:    os << "indefinite";     break;
      case Missing:       os << "missing";        break;
   }
   return os;
}



// boolean.hxx ------------------------------------------------------------

#ifndef BOOL_LIB
#define BOOL_LIB 0

enum BOOL { FALSE, TRUE };

#endif


// end ----------------------------------------------------------------------